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PREFACE 
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The  work  described  is  part  of  an  investigation  of  the  propagation  ol  elastic  waves  in  soil 
materials  near  the  freezing  point  of  water  sponsored  by  ihe  Advanced  Research  Projects 
Agency  (ARPA  Order  152!>).  This  report  deals  with  efforts  to  measure  cotnpiosnioual 
and  shear  velocities,  and  their  respective  attenuations,  from  the  modes  of  free  vibration  of 
an  elastic  sphere.  An  attempt  has  tieen  made  to  provide  a  reasonably  complete  and  self- 


contained  account  of  the  appropriate  theory  (much  of  which  is  tutorial),  a  useful  description 
of  the  experimental  and  data-reducing  processes,  and  some  graphical  aids  to  allow  manual 
data  reduction.  Computer  codes  are  not  listed  but  may  he  obtained  (in  Fortran  11)  from  the 
author. 
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ON  THE  DETERMINATION  OP  ELA8T1C  AND  ANELA8TIC 
PROPERTIES  OP  ISOTROPIC  SPHERES 

by 

Martin  L.  Smith 


INTRODUCTION 

An  elastic  sphere  whose  properties  vary,  at  most,  as  a  function  of  radius  alone  is  the  only 
bounded  three-dimensional  body  whose  entire  spectrum  of  free  vibrations  can  be  determined  exactly. 
(A  quantity  is,  here,  considered  to  be  "determined  exactly"  if  it  can  be  computed  as  the  root  of 
an  expression  which  is  composed  of  a  finite  number  of  conventional  transcendental  functions  or 
which  can  be  evaluated  by  a  one-dimensional  numerical  integration.)  Roughly  speaking,  this  is  a 
consequence  of  the  fact  that  in  any  of  the  three  elementary  coordinate  systems  (Cartesian,  cylindrical 
and  spherical  coordinates)  a  sphere  of  the  above  composition  is  the  only  finite  body  which  is  totally 
symmetric  under  two  of  the  three  coordinates. 

Because  the  solutions  so  obtained  are  not  restricted  to  extremely  high  or  low  frequencies  or  to 
limiting  geometries,  as  in  the  case  of  cylinders  and  plates,  we  may  reasonably  hope  to  evolve  an 
interferometric  technique  for  inferring  compressional  and  shear  velocities  from  observed  resonance 
spectra.  This  has,  in  fact,  been  done. 

The  earliest  such  work  is  that  of  Birch  (personal  communication)  which  antedates  World  War  II, 
and  is  unpublished.  The  earliest  published  work  is  that  of  Fraser  and  LeCraw  (1964).  This  pioneer 
work  was  done  on  small  (19-24  mm  radius)  spheres  of  yttrium  gallium  garnet  and  yttrium  aluminum 
garnet.  Fraser  (1968,  1970)  made  extensive  use  of  the  technique  to  study  the  properties  of  vitreous 
silica.  The  techniques  presented  by  Fraser  and  LeCraw  (1964)  were  also  applied  by  Anderson  and 
Soga  (1966)  to  very  small  (0.5  mm)  spheres  of  polycrystalline  MgO,  by  Soga  and  Anderson  (1967)  to 
two  tektites,  moldanite  and  indochinite,  and  by  Anderson  et  al.  (1970)  to  very  small  glass  spheres 
(about  0.5  mm)  taken  from  a  sample  of  lunar  fines. 

The  use  of  spherical  interferometry  to  infer  elastic  properties  requires  that  we  be  able  to 
"invert”  a  known  geometry  and  a  set  of  observed  resonant  frequencies  into,  in  the  case  of  a 
homogeneous  sphere,  a  compressional  and  a  shear  velocity.  Previous  workers  have  used  graphical 
techniques  based  upon  forward  eigenfrequency  calculations  by  Cole  and  Fraser  (unpublished).  Al¬ 
though  such  methods  do,  in  fact,  yield  high  accuracy  for  geometrically  precise,  high-Q  materials 
such  as  vitreous  silica  (Fraser  1968),  we  believe  that  the  development  of  statistical,  relatively 
objective  methods  was  appropriate  for  the  study  of  much  less  ideal  materials,  such  as  frozen  soils. 

In  later  sections  we  describe,  and  use,  a  numerical  algorithm  for  inverting  arbitrary  sets  of  data  and 
describe  a  rough  method  of  estimating  uncertainties  in  the  result. 

A  second  obstacle  we  have  encountered  in  the  use  of  low-Q  (Q  <  500)  materials  is  the  difficulty 
of  identifying  modes,  since  many  resonances  are  lost  in  the  "skirts”  of  others.  This  difficulty  is 
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less  pronounced,  and  uay  often  be  absent,  with  high-Q  materials.  In  relief  of  this,  we  have 
developed  the  appropriate  theory  of  the  forced  response  of  an  elastic  sphere  and  utilized  it  to  esti- 
mate  relative  modal  amplitudes.  These  results  have  been  extremely  useful  in  interpreting  experi- 
mental  results. 

Wo  have  also  extended  both  the  forward  and  inverse  algorithms  to  deal  with  arbitrarily  layered 
spheres  (which  still,  of  course,  retain  spherical  symmetry),  in  order  to  admit  of  jacketed  samples. 
Computer  storage  limitations  presently  prohibit  multilayer  inverse  calculations  and  we  do  not 
present  here  numerical  or  experimental  results  on  such  samples. 

In  pursuit  of  various  of  the  above,  we  develop  expressions  for  the  kinetic  and  various  potential 
energy  densities  associated  with  modes  of  free  vibration  and  show  how  these  are  useful  in  a)  esti¬ 
mating  the  effect  of  compositional  perturbations  on  the  eigenspectrum  and  b)  partitioning  a  given 
mode’s  lumped  Q  into  compressional  and  shear  Q’s.  Coquin  (1964b),  in  an  application  of  his 
earlier  technique  (Coquin  1964a),  has  also  solved  the  latter  problem  for  a  homogeneous  sphere  by 
using  its  characteristic  equation,  His  results,  while  obtained  by  a  different  technique,  are 
fundamentally  identical  with  ours,  and  could  be  extended  numerically  to  handle  the  multilayered 
version  of  either  a)  or  b)  above, 

Many  of  the  techniques  and  results  we  have  utilized  are  drawn  from  seismological  literature 
concerned  with  interpretation  of  the  earth's  free  oscillations,  Particularly  lucid  theoretical  dis¬ 
cussions  are  available  in  Backus  and  Gilbert  (1967)  and  Dahlen  (1968).  The  methods  detailed  by 
Backus  (1967)  are  extremely  helpful  in  problems  of  this  sort  and  we  have  made  extensive  use  of 
them.  One  of  the  most  useful  discussions  of  the  behavior  of  vibrating  systems  is  presented  in 
Rayleigh’s  Theory  of  Sound  (Rayleigh  1945),  first  published  in  1877. 


ELASTIC  DISPLACEMENT  SOLUTIONS  IN  SPHERICAL  COORDINATES 

We  consider  a  volume  of  space  filled  with  an  isotropic,  homogeneous,  linearly  elastic  medium 
having  Lamd  constants  A  and  p,  and  density  p.  Vie  assume  the  medium  to  be  free  of  gravitation  and 
other  body  forces,  but  allow  the  existence  of  one  of  more  surfaces  across  which  tractions  may  be 
applied. 

Let  u  be  the  displacement  field  specifying  the  motion  of  each  particle  from  its  unique  rest 
position.  We  assume  n  to  be  a  first  order  infinitesimal  and  do  not,  in  the  absence  of  zeroth  order 
fields,  have  to  distinguish  between  Eulerian  and  Lagrangian  coordinate  systems.  Let  T  be  the 
elastic  stress  tensor.  If  we  assume  that  u  =  0  corresponds  to  the  unstressed  state  of  the  medium, 
T  is  given  by  (Fung  1965): 


T  =  A(V»  u H  +  p(Vu  +  uV),  (1) 

where  £  is  the  identity  tensor.  Vu  is  the  gradient  tensor  of  u,  and  uvis  its  transpose. 

The  conservation  of  linear  momentum  leads  immediately  to  the  equation  of  motion, 

pdf  u  =  V  •  T.  (2) 

Equation  1  and  some  standard  vector  calculus  identities  convert  eq  2  to 
pdfu  =  (A  +  2p)V(V  •  u)  -  pV  x  V  x  u. 


(3) 
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We  choose  to  represent  u  by  (Backus  1967) 

h  =70  +  VYV  -Tx  VtW  (4) 

where  U,  V,  and  W  are  scalar  fields  and  vt  is  defined  by 

+  <^(sin  6)~ld<j).  (5) 

Tis  a  unit  vector  directed  away  from  the  origin,  d  and  </>  are  the  colatitude  and  longitude,  and  6  and 
<£are  their  respective  unit  vectors.  Vj  is  the  gradient  operator  on  the  surface  of  a  sphere  of  unit 
radius.  It  is  related  to  the  three-dimensional  gradient  operator  by 

V  =  Tdt  -  r-lVt.  (6) 

After  seme  algebra,  we  can  show  (Backus,  1967)  that 
V(V  •  tt)  =Tar  |(<?r  +  Z-)u  +  rlV?p}  + 

+  +  ^)U  +  r*VlV}'  (7) 


and 

V  x  v  x  u  =  T\r~zdtixv\V)  -  r_2v2U i  +  Vilr-1<?rl/  -  r~ld2(rV) I  + 

+  T  x  ir-2vfw  +  r-1<?f  Wl.  (8) 

We  insert  eq  4,  7  and  8  into  eq  3.  We  now  appeal  to  the  uniqueness  of  the  representation  4  (Backus 
1967)  to  yield  the  three  coupled  partial  differential  equations 


Pd2tU  =  (A  +  2(i)dI  |(<?r  +  i)  U  +  rlv^}  - 

-  I  dtij^V)  -  •  (9a) 

r2 

pdfv  =  (A  +  2p)  |^r-1<?r  +  ~)  U  +  ^lV}  “ 

-  ^  \drU  -  d2itV)\,  (9b) 

r  r 


and 


pd2W  =  t  )dr2(rW)  +  rVi^l- 


(9c) 
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We  note  that  both  V  and  W  may  be  augmented  by  any  constant  without  affecting  n  (see  eq  4).  There¬ 
fore,  we  may  expect  these  two  scalar  fields  to  be  determined  only  to  within  an  additive  constant. 

We  also  note  that  eq  9a  and  9b  both  involve  V  and  V  but  not  W,  while  eq  9c  involves  W  but  neither 
1/  nor  V, 

The  manner  in  which  one  elects  to  solve  eq  9a,  9b  and  9c  (plus  whatever  boundary  conditions 
appertain)  depends  upon  the  intended  application  of  the  solution.  For  our  purposes,  we  wish  to  find 
a  set  of  linearly  independent  vector  fields,  each  of  which  satisfies  eq  3.  If  the  set  is  complete, 
all  possible  solutions  to  eq  3  may  then  be  expressed  as  some  linear  combination  of  the  members  of 
the  set,  the  coefficients  used  in  the  expansion  being  determined  by  boundary  and  initial  conditions, 

In  pursuit  of  this,  we  Fourier  transform  eq  9a-9c,  going  from  time  t  to  angular  frequency  w.  We 
do  not  introduce  a  distinct  symool  for  Fourier  transforms  since  it  will  be  clear  from  the  context 
whether  a  symbol  refers  to  the  transformed  or  untransformed  variable.  The  result  of  Fourier  trans¬ 
formation  is  to  replace  d2  by  -a2. 

To  transform  the  i  isultant  tri '  of  equations  from  partial  to  ordinary  differential  equations  we 
introduce  a  surface  spherical  harmonic  expansion  of  U,  V  and  W.  For  I  >0  and  -l  <  m  <  l,  we 
define  (HiU  1953) 

Y?  %  <f>)  =  (-1V"  f^L tJ  1-Z  jmj|1|“  P®  (cos  0)e'®<*  (10) 

I  4 n  U  +  \m\)\  l 

where  P®  is  the  Associated  Legendre  Function  given  by 

Pj"(x)  =  (Li_»8)(|m|/8)  ji+|a|  0,2  _  1}It  (n) 

2*1! 

If  Sj  is  the  surface  of  a  sphere  of  unit  radius  centered  on  the  origin,  the  F®  are  orthonormal  in  the 
sense 


/  F®  (0.  *)  Fft  (0,  4)  sin  0d0d<^  =  S]XSmil  (12) 

S1 

where  the  bar  indicates  complex  conjugation  and  <5^  is  the  Kronecker  delta. 

The  F®  (0,  <f))  form  a  complete  set  and,  if  we  assume  that  U,  V  and  W  are  sufficiently  regular, 
we  may  expand  them  as 

OO  l 


U  (r,  0,  <j),  w) 

V  (r,  0,  4>,  «) 


1  1  (/™(r,  <u)F®  (0,  0), 

1=0  m=— 1 


I  i  F®(r,  «)F®  (0.  $), 
1-1  m=-l 


(13a) 


(13b) 


and 


W  (r,  0,  0,  <o) 


1  i  W?  (r,  w)  Ff  (0,  <f>). 

1=1  m=— 1 


(13c) 


The  i  =  0  terms  have  been  omitted  from  eq  13b  and  13c  since  inspection  of  eq  4  reveals  that  these 
terms  do  not  contribute  to  the  displacement  field. 
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We  now  insert  eq  13a- 13c  into  the  transformed  equations  9a-9c,  and  make  use  of  the  relation 

*=  -  ni  +  i)y“.  (14) 

The  resultant  expressions  are  then  multiplied  by  a  particular  Yf  and  integrated  over  the  surface  of 
a  sphere  of  radius  r.  Application  ol  eq  12  leads  immediately  to  the  set  of  coupled  ordinary  differential 
equations: 

-  pu‘ii f  .  (A  *  ws,  {(a,  *  J )of  -  17}  * 

+  JL  ||(I  +  l)a,(rV”)  -  HI  +  l)Uf\  (15a) 

2  r 
r 

-  .  <i  *  ad  {('■»,  *  ^)"r  - 


-  ii  |5rl/“  -  8*(rVpl. 


(15b) 


and 

_  (rW™)  -  IU  W”!  .  (15c) 

We  note  that  none  of  eq  15a-15c  explicitly  involves  m. 

The  set  15a-15c  can  be  solved  by  any  of  several  standard  techniques.  The  method  used  here 
is  detailed  in  Appendix  A.  The  results  are 


'd,j\(kr) 

lx(kr) 


Vlttr) 

yx  ( kr ) 

r-1<Jr[r;,(yr)]  — — 


HI  +  1) 


yi(yr) 


‘^Ir/j  (yr)] 


for  l  >  0; 


U 


0 

0 


\dtj0  (kr) 


d/oOcrH 


(16) 

(17a) 


and 


y°o  -  0; 

H'f  =  Ir/jlyr) 


for  1  >  0; 


(17b) 

(18a) 


(18b) 
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where 


(19) 


(20) 


V  and  Vn  ure,  respectively,  the  velocities  of  compressions!  and  shear  waves  in  the  medium.  The 
functions  /,(x)  and  y((x)  are  spherical  Bessel  functions  defined  by 


where  J,  and  V,  are  conventional  cylindrical  Bessel  functions. 

We  note  that  the  functional  form  of  the  solutions  16*18  is  independent  of  m.  The  degeneracy  of 
this  index,  which  will  persist  to  the  case  of  a  layered  sphere,  is  a  direct  result  of  spherical  symmetry. 
Any  linear  combination  of  surface  spherical  harmonics  of  order  1, 


tf,  =  k  a™  K™  (23) 

m=-l 

is  itself  a  surface  spherical  harmonic  of  order  i  and  satisfies  eq  14. 


FREE  OSCILLATIONS  OF  A  LAYERED  ELASTIC  SPHERE 

We  consider  a  sphere  divided  into  N  concentric  spherical  shells.  We  number  the  shells  from 
the  center  outwards  and  let  Tj  be  the  outermost  radius  of  the  ltfl  shell.  Then  rN  is  the  radius  of  the 
sphere.  Let  r0  equal  zero.  We  suppose  that  the  Ith  shell  is  composed  of  an  elastic,  homogeneous, 
isotropic  medium  of  density  pt  and  having  Lamd  parameters  Aj  and  gj.  These  parameters  define  the 
congressional  velocity  Ppl  and  the  shear  velocity  Vgl. 

We  assume  that  the  surface  r  =  rN  is  free  from  all  tractions,  and  that  no  body  forces  (such  us 
gravitation)  are  present.  We  wish  to  know  for  what  angular  frequencies  w  there  exists  a  non-trivial 
displacement  o(r)  e,cul,  such  that  the  surface  is  free  of  traction  and  all  internal  boundary  conditions 
(discussed  below)  are  fulfilled.  We  shall  label  such  angular  frequencies  eigenfrequencles  and  their 
associated  displacements  eigenlunctions.  We  refer  to  such  traction-free  motions  as  free  oscillations. 
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Lot  u(,,U)  .>*'  l»o  tlm  displacement  field  in  the  Ith  layer.  From  tin*  displacement.  we  can  compote 
a  stress  tensor  T*1*  by  et|  1.  Equation  2,  namely 


V  •  T 


(2) 


must  hold  everywhere  in  the  meditnn,  since  it  excesses  only  the  conservation  of  linear  nionientoin  and 
is  not  dependent  iijkiii  such  assumptions  as  isotropy,  homogeneity,  etc.  In  particular,  eq  2  is  valid 
in  a  small  "Gaussian  pillbox"  which  encompasses  a  portion  of  the  surface  r  r, .  the  boundary 
between  the  ilh  and  (I  *  l),,,  shells.  Let  i  denote  the  volume  enclosed  between  the  radii  r, 
and  r,  ♦  Sr  and  by  some  ratine  of  the  coordinates  0  and  6.  Then 


Fa 


(  pd~  yde  ■  /  T dr. 


(24) 


If  1  is  the  surface  of  o.  Gauss's  theorem  leads  to 
/  pdfiide  -  J  T  •  n 'do 

V  1 

where  a  is  the  unit  outward  normal  on  1 .  If  we  let  Fa  go  to  zero,  then  the  volume  or  c  also  vanishes 
and,  if  remains  bounded,  the  left-hand  integral  goes  to  zero.  The  right-hand  side  becomes 


JIT'**0  -  T<l“,)  I  •  Ydti  0  at  r  -  r, .  (25) 


Since  this  result  ts  independent  of  the  details  of  the  shape  of  '  we  conclude  that  the  quantity  T_  •  r 
must  be  everywhere  continuous,  and  in  particular  across  boundaries. 

To  condition  25  we  add  one  expressing  our  intuition  of  the  behavior  of  elastic  materials.  If 
both  the  i,ft  and  (/  ♦  l)1'1  shells  are  solid,  we  require  that  the  displacement  u  be  continuous.  We 
refer  to  interfaces  at  which  this  is  true  as  being  "welded."  If,  however,  one  or  both  shells  are 
fluid  (that  in,  p  0),  we  require  only  the  radial  component  of  displacement  to  bo  continuous.  In  the 
latter  case,  we  allow  the  boundary  to  slip  laterally  but  in  neither  case  do  we  allow  holes  to  open  or 
matter  to  interpenetrate  itself. 

The  quantity  T  •  7is  a  vector  and  is  the  traction  (force)  acting  on  a  surface  ncnnal  to 
r7  in  a  fashion  identical  to  eq  4,  we  may  represent  I*  as 

T  .  7*  -  ?P  *  VtQ  -  7  •  V,R.  (26) 


where  P,  Q  and  R  are  scalar  fields.  Equation  25  states  that  P,  Q  and  R  are  continuous  across  an 
interface.  The  stress-displacement  relations,  eq  1,  enable  ns  to  relate  P,  Q  and  R  to  U,  V  and  K, 
the  scalar  representatives  for  a.  These  relations,  which  are  derived  in  Appendix  B,  arc 


P 


(A  ♦  2p)  r)f  V  +  -~U 


(27a) 


(27b) 
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end 


*  •  «'*,( *)• 


(27c) 


The  import  of  the  boundary  conditions,  then,  is  that  P,  Q,  R  and  U  aro  continuous  everywhere 
and  V  and  W  are  continuous  in  solid  domains. 


We  now  have  six  scalar  fields  to  contend  with.  Hov  ever,  an  examination  of  eq  16,  17,  18  and 
27  indicates  that  we  can  group  them  into  two  sets,  one  consisting  of  U,  V,  P  and  Q,  and  ono 
consisting  of  W  and  R.  These  two  sets  are  completely  independent;  they  do  not  interact  in  any 
way.  We  may.  without  loss  of  generality,  treat  them  separately. 

The  set  (U ,  V,  P,  Q)  is  the  set  of  spheroidal  variables.  The  displacements  described  by  this 
set  give  rise  to  both  compressional  and  shear  strain,  and  are  interference  products  of  compressional 
and  vertically  polarized  shear  waves.  Rayleigh  surface  waves  are  one  product  of  this  class  of 
motion. 

The  set  (W,  R)  is  the  set  of  toroidal  variables.  The  displacements  described  by  this  set  are 
orthogonal  tor  and  produce  only  shear  strain.  They  are  the  interference  products  of  horizontally 
polarized  shear  waves.  Love  surface  waves  are  associated  with  this  set. 

From  this  point  forward  we  will  discuss  only  spheroidal  types  of  motion.  A  similar  development 
for  toroidal  oscillations  can  be  easily  formulated  since  the  toroidal  problem  is  substantially  simpler. 

For  any  specified  angular  frequency  of  oscillation  to,  the  set  (17.  V,  P,  (?)  can  be  expanded  in 
terms  of  spherical  harmonics  as 


U  (r,  0,  (o)  = 

<x 

V 

1=0 

— II 

B 

(7®  (r,  a)  K®  ( 0 ,  <J) 

(13a) 

V  (r,  0,  <£.  to)  m 

rx 

V 

1=1 

t 

m=-l 

F®  (r,  to)  V®  (0,  6) 

(13b) 

P  (r,  0,  <f>,  to)  = 

IX 

V 

UO 

i 

m--I 

Pf(r,  to)  V®  (0,  6) 

(28a) 

Q(r,  0,  <4,  to)  = 

ix 

V 

1  1 

m  -I 

Qf  (r,  <o)  F®(0.  tf). 

(28b) 

The  U ®  and  F®  ate  ven,  in  terms  of  a  set  of  coefficients,  by  the  analytic  solutions  16  and  17. 
P ®  and  Q®  are  then  given  by  eq  27.  Our  problem  is  to  determine  those  frequencies,  to,  for  which 
we  can  construct  U,  V,  P  and  Q  for  each  layer,  such  that  all  boundary  and  interface  conditions  are 
met. 

Because  the  V®  are  au  orthogonal  set.  we  may  consider  the  above  problem  separately  for  each 
F®.  That  is,  given 


U(r,  0,  4,  to)  =  If®  (r,  o,)  Yf(0,  <4),  (29) 

etc.,  for  what  angular  frequencies  o»  can  we  satisfy  all  internal  aud  external  boundary  conditions? 
Since  m,  as  discussed  earlier,  is  u  degenerate  index,  we  can  simplify  eq  29  to 
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uo.o.t.w)  -  t/,  (r.  oj)  M,  (0.  $), 

(30a) 

V(r.0,4.*>)  -  y,(r.M)N|(0.4). 

(30b) 

P(r,  0,  *.„.)  P,(r. «)«,(».  4). 

(30c) 

Q{t.  0.  tf.  o,)  g,(r.  o>)  //,  si). 

(30d) 

whore  1/ 1  is  given  by  eq  23  and  is  some  spherical  harmonic  of  degree  /.  *  he  details  of  //|  are  of 

no  particular  interest  at  present  since  the  clgonfreqnencies  ;uid  the*  form  of  Ut . <?,  remain  fixed  no 

matter  what  H,  we  choose. 

Otir  problem  now  is  to  find  the  eigenfreqitencies  associated  with  spherical  harmonics  of  degree 
/.  If  we  wish  to  know  all  e'genfrequencies.  we  must  repeat  this  procedure  for  each  of  /  0.1. 

•) 

s»*  •  •  • 

We  will  devise  a  constructive  algorithm  which  will  enable  ns,  Tor  a  given  angular  frequency  <.j. 
to  construct  u  solution  from  the  center  outward  which  meets  all  boundary  conditions  save  one.  If 
the  last  condition  is  met.  <.»  is  an  eigenfrequency.  The  method  we  present  here  is  one  application  of 
the  use  of  propagator  matrices.  The  general  method  is  discussed  by  Gilbert  and  Backus  (1966). 
Seismological  applications  are  numerous  (see,  for  example.  Harkrider  196-1.  Ben-Menahem  1904b  and 
Anderson  1966). 

We  consider  eq  16.  which  expresses  U,  and  V,  as  a  linear  combination  of  four  independent  func¬ 
tions  of  radius.  Equation  27  allows  us  to  extend  eq  16  to 


Vi 

P{ 

<?{ 

V  J 


B. 


(31) 


where  /  »  1 . N  designates  the  layer  to  which  the  solution  is  appropriate.  //{ is  a  4  \  4  matrix 

constructed  from 


.  o  \  j  i,  2A  .  .  f(i  +  1)  . 

)3j  .  (A  ♦  2fj)dr/t,j  +  —  A|j  -  A - - - Agj, 


(32a) 


and 

h4j  -  ♦  rryr'/tgj)!  (32b) 

and  the  first  two  rows  of  H\  are  taken  from  eq  16.  The  set  M{ . D{l  are  constants.  We  omit  the 

0  and  ^  terms  for  convenience. 

We  rewrite  eq  31  as 

«'(r)  -  HJ(r)  •  C‘  (33) 

where  wo  have  omitted  the  subscript  J.  The  vector  S*(r)  includes  both  the  stress  and  displacement 
terms. 
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We  will  now  prrx'ccd  to  construct  a  solution  satisfying  all  internal  boundary  conditions.  In 
region  I  which  includes  the  point  r  0  we  can  a  prmn  eliminate  those  solutions  winch  go  as 
yt  (Ar)  and  y((yr).  Therefore,  C1  has  the  fonn 


where  e*.  and  e7,  are  Euclidean  four-dimensional  unit  coordinate  vectors  directed  along  the  first  and 
second  coordinate  axes.  In  the  first  region,  then,  we  have 

S*(r)  =  ff'(r)  •  U'e*  t  fl'e'l.  (36) 

lu  the  second  region,  we  have 

S~ (r)  =  Hs(r)  •  C2.  (36) 

Assuniing  both  regions  to  be  solid,  the  boundary  conditions  require  that 

S2(r,)  =  S  l(r,).  07) 


or 


»2(r j)  •  C2  // 1  (r,)  •  M'e*,  +  Bl e^l.  (38) 

Because  the  solutions  composing  the  columns  of  Rare  linearly  independent,  matrix  theory  guarantees 
lhat  H  is  nonsingular.  Therefore,  we  may  express  C2  as 

C2  I /(2(r,)|  1  •  H 1  (r j )  •  1/t‘e*  +  fl'e*!.  (39) 

An  alternative  form  for  oq  39  is 

C2  A'{~  *  »'<"  (40) 

whore 

£2  =  i Hs(rl)|-‘  •  U\r j)  •  e*  (41) 

£a  l//2(r,)r'  •  }l 1  ( r j )  .  (42) 

Equations  40,  41  and  42  suffice  to  specify  S2(r).  By  a  similar  procedure,  we  can  extend  the 
solution  from  the  ilh  to  the  (i  *  l),ft  shell.  The  appropriate  relations  are 

Sifl(r)  .  l[hHr)  •  C*'1  (43) 

CU1  A1?"  i 


(44) 
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£nt  =  |//' 1 1  (r )  1“ 1  •  fl‘(r)  •  (i:,) 

«",f  .  |//itf(r)l  *  •  II1  (r)  ■  £' .  (4b) 

In  fact,  either  of  * 1  or  1  *  can  be  computed  alone  by  starting  with  eq  41  or  42  and  simply 
applying  eq  45  or  46  as  many  times  as  is  necessary. 

We  suppose  now  that,  we  have  computed  (N  iutd  £N  wliore  N  is  the  number  of  shells.  The 
solution  in  this  shell  is  given  by 

S N(r)  HN(r)  •  \AHN  •  B,CW  1.  (17) 

We  note  that  both  the  solution  associated  witli  li.e.,  (r)  when  B^  ()|  and  the  solution 
associated  with  £N  separately  satisfy  all  of  the  internal  boundary  conditions.  If  w  is  an  eigen- 
lrequency.  we  will  be  able  to  find  some  .4*  ami  B  ’,  not  both  zero,  such  that  the  last  two 
components  of  S N  (r^)  are  zero. 

A  straiglitforward  way  to  do  this  is  to  evaluate 

(a)S*  =  JlN  (rN)  •  £N,  ^ 

and 

<b>Sw  -  HN(rN)  ■  £N.  (48» 

If  S*  (rN),  i.e.  P,  must  vanish,  A 1  and  B1  must  satisfy 

AHaS")  =  -  Bl(bS").  (49) 

We  may,  without  loss  of  generality,  impose  a  normalization  such  as 

U41)2  f  (B1)2  =  1.  ^,0) 

which,  with  eq  49,  allows  us  to  compute  -41  and  B1.  Using  A 1  and  B1,  we  then  compute  (rN),  i.e. 
Q,  and  examine  it  to  see  if  it  vanishes.  If  it  does,  u>  is  an  eigenfrequency;  it  it  does  not,  m  is  not 
an  eigenfrequenev. 

If  oj  is  an  eigenfrequency  then  we  may  use  the  values  of  -41  and  B1  in  eq  4.5  and  44  to  compute 
the  eigenfunction  S *(r)  at  any  radius  in  the  system.  We  recall  that  S1  (r)  must  be  multiplied  by  some 
spherical  harmonic  of  degree  /  and  the  factor  to  obtain  the  full  solution, 

S{  (r,  0,  </;.  0  S|(r) //,((?,  <f>)  eiwt 

where  we  have  reinstated  the  subscript  /.  We  emphasize  again  that  the  precise  nature  o.  //,  depends 
upon  other  considerations  and  is  not  relevant  here. 

This  method  requires  modification  when  either  l  0,  or  /  /  0  but  one  or  more  shells  have  a 
vanishing  shear  modulus.  We  will  briefly  outline  the  fonn  these  modifications  take. 

When  /  =  0,  B1  vanishes  identically  and  only  one  solution  is  propagated.  The  matrix  j/j  is 
collapsed  to  a  2  ’  x  2  matrix  by  eliminating  those  solutions  with  arguments  yr.  V  and  Q  both  vanish 
and  the  only  condition  at  r  =  rN  is  that  P  vanish.  A]  becomes  merely  a  scale  factor  and  may  be 
taken  as  equal  to  unity. 
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In  a  fluid  region,  Q  vanishes  identically.  Across  a  solid /fluid  or  fluid/fluid  interface  V  may  he 
discontinuous  and  Q  is  continuous  and  zero.  If  we  are  propagating  upward  through  a  solid  and 
encouuter  a  fluid,  we  must  combine  the  £  and  £  solutions  to  yield  Q  0  at  the  interface.  This  re¬ 
quirement  determines  A1  and  B1  and,  therefore,  &  for  all  shells  up  to,  and  including,  the  fluid 
region. 

Consequently,  in  a  fluid  region  we  have  only  one  solution.  In  crossing  from  a  fluid  to  a  solid 
region,  we  must  "start”  a  new  solution  having  some  non-zero  V,  but  for  which  U ,  P  and  Q  vanish, 

We  may  always  find  some  £  which  yields  this  result. 

For  a  given  1,  we  shall  arrant, e  the  frequencies  of  free  oscillation  in  ascending  order,  as  qo^, 
jWj,  .jOj | ,  ...  We  shall  designate  the  displacements,  as  a  function  of  r,  as  0u,,  jUj,  ...  and  the 

four -vector  of  displacement  and  stress  by  0Sj .  jSj . etc.  The  "lowest”  mode,  for  a  given  1  is 

referred  to  as  the  "fundamental"  mode  and  the  remainder  as  "overtones." 

Each  mode  of  order  l  represents,  in  fact,  a  space  of  modes  of  dimension  21+1  (the  number  of 
possible  values  for  m  between  -1  and  +1).  Each  member  of  this  space  has  the  same  eigenvalue  and 
depends  functionally  upon  radius  in  the  same  way,  but  is  associated  with  a  different  surface 
spherical  harmonic  of  order  1.  The  entire  set  of  modes  for  a  given  layered  sphere  consists  of  the 
doubly  infinite  sets  (0  <  n  <  <*>  and  0  <  1  <  <~)  of  21+1  degenerate  spheroidal  and  toroidal  modes. 

This  set  is  complete  (Backus  1967)  and  therefore  any  elastic  motion,  consistent  with  the  fore¬ 
going  boundary  conditions,  which  the  system  may  undergo  can  be  described  as  a  sum,  in  the  time 
domain,  of  members  of  this  set.  A  discussion  of  the  asymptotic  relations  between  modes  and 
"waves"  is  given  by  Ben-Menahem  (1964).  As  an  example  of  this  we  note  that 

Y\lfi 

as  can  be  seen  from  eq  10  and  11.  Y{  describes  motion  which  is  closely  confined  to  the  equator 
(0  n/2)  and  which  behaves  as  waves  traveling  circumferentially  about  the  equatorial  zone.  The 

change  of  ptiase,  per  unit  of  distance  traveled  in  the  ^direction,  is  1/a  and  is  therefore  the  mode’s 
surface  wave  number.  For  large  values  of  cu,  we  may,  rough'v,  expect  that  the  quantity 

U*  =  a  —  (52) 

ell 

represents  a  group  velocity  and 

C*  =  —  (53) 

I 

represents  a  phase  velocity, 

ENERGY  DENSITIES  AND  RAYLEIGH’S  PRINCIPLE 

The  nth  overtone  of  the  spheroidal  mode  of  degree  I  has  displacements  nu,,  given  by 

n*i  =  •?*„(/, (r)  «,((>,<£)  +  nFj(r)  V,H,(0.<£) 


,  nf  hi  +  1  (sin  eii4> 


4ct(2I)! 


(54) 
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where  H j  is  some  surface  harmonic  of  decree  /.  The  left  subscript,  n,  donolos  tin*  overtime  to  which 
U  and  V  are  appropriate. 

The  total  kinetic  energy  of  motion,  averaged  over  one  cycle,  is  given  by 


E„  =  I  -2 


K  =  o  nwl  /  PW  •  Wdv‘ 


(55) 


Using  the  results  of  Backus  (1967),  we  can  show 


4  0 


(56) 


when  W[  is  normalized  such  that 


/  Wj/fj  sin  0d0d<, 5  =  1 


(57) 


(see  eq  12). 

The  total  elastic  energy  is  given  by 


1 


Ea  =  i/tr(T  •  e)  dV 


(58) 


where  lr  indicates  the  trace  operator,  T  is,  as  before,  the  elastic  stress  tensor,  and  e  is  the 
(infinitesimal)  strain  tensor  defined  by 


e  =  -  IVu  +  uVl. 
2 


(59) 


After  some  algebra,  we  can  show 


Ea  =  i  f  j  At  r  r?r  U  +  2  V  -  1(1  +  1)V]2  +  p[2(rry/)2  f  1(1  +  1  )(U  +  rd/  -  V)2 


[2 U  1(1  +  1)V]2  +  (l  -  !)/(/  +  1)U  *  2)V2]}dr 


(60) 


where,  for  simplicity,  we  have  dropped  the  subscripts  n  and  /. 


When  ntij  is  an  eigenfunction  of  the  system  we  are  interested  in,  and  najJ"  is  its  associated 


eigenvalue,  we  have 


EK  -  Ea‘ 


(61) 


If  we  define 


K  =  [(/2  c  1(1  +  DV2]  r2, 


(62) 
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L  \r<\U  t  2 U  -  1(1  i  1)1/!",  and  (63) 

M  2(r,\U)2  i  1(1  f  1)([/  ,  r<\V  l/)"  ( 

f  [2U  -  1(1  t  1)F|2  ,  (1  1)1(1  f  1)(1  ,  2)V2  (64) 

we  may  rewrite  eq  61  as 

rN 

/  (AL  f  pM)dr 

o  0 

„^r  = - •  (65) 

rN 

/  pK  dr 
o 

As  first  pointed  out  by  Rayleigh  (1945)  and  discussed  by  Backus  and  Gilbert  (1967),  the 

g 

functional  n«q  (uU|),  defined  by  eq  62-65,  is  stationary  with  respect  to  small  variations  in  U|  if 
nu,  is  an  eigenfunction  associated  witlt  As  a  consequence  of  this,  we  may  use  eq  65  to  relate, 
to  first  order,  small  changes  in  A,  p  and  p  to  small  changes  in  t|i up,  The  result  is 


rN 

f(SAL  ,  SpM  Sp„rK)dr 
0 

rN 

f  pKdr 
o 


wltere,  clearly, 


) 


*(„“.>■ 


(66) 


The  ‘ 'partial  deiivatives"  of  nm|  with  respect  to  A,  p  and  p  may  be  converted  to  those  with  respect 
to  V  ,  V'j,  and  p  by  application  of 


V2 

\  tin  /  „  „  s 

(67a) 

WVp'Vs 

f(U\  v2  -  2V2 

U/vt.v  " 

(67b) 

p’  s 

(ir)  ■  2pV° 

(67c) 

(— \  -  -4 pV 

UJvD.p 

(67d) 
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(67e) 


(67f) 


EFFECTS  OF  A  SLIGHT  ANELASTICITY 

We  now  suppose  that  the  perturbations  in  A  and  p  aie  purel”  imaginary: 


u  A  /(jjA1 

(68a) 

op  iio/i' 

(68b) 

where  A1  and  p1  are  real.  This  is  precisely  equivalent  to  replacing  A  by  A  +  AV?t  and  p  by 
p  4  H‘At  in  eq  3,  Then 


*(nwl 


rN 

i  ntU|  f  (A'L  4  [i‘M)(lr 
o 

rN 

/  i>Kdr 

o 


(69) 


to  first  order  in  A1  and  p\  Alternatively, 


a(n"P  =  J  n<tJl  nDl 


which  serves  to  define  Dj.  Therefore 


g  ,D1 

and  the  lumped  quality  factor  nQ,  is  given  by 


Q,  — 

"  *  nDl 


We  define 


Q 


-t 


p  A  4  2  p 


(A'  +  2p') 


and 


°v 


(70) 


(71) 


(72) 


(73) 


(74) 


16 


ELASTIC  and  an  elastic  properties  of  isotropic  spheres 


These  are,  respectively,  the  inverse  Q' s  of  pure  congressional  and  pure  shear  waves  in  a  homogene¬ 
ous,  slightly  lossy  medium.  From  eq  76  and  74  it  follows  that 


Combining  eq  69,  70,  72,  75  and  76,  we  have 


n 


-1 


X 0 


IQ"?  (A  t  2fi)L 
1> 


Q'J  (tiM  -  2gL)Jdr 

v  S 


fN 

f  pK  dr 

0 


(75) 

(76) 


or 


n 


<?;1 


rN 

/  [Q'vl  (A  f  2p)L  f  <?*'  f*«  "  2L)] dr 
0  l» 


2£k 


(77) 


RESULTS  FOR  A  LAYERED  SYSTEM 

The  above  expressions  are  all  valid  for  a  spherically  symmetric  body  whose  properties  A,  p 
and  p  are  arbitrary  (but  reasonable)  functions  of  position.  For  the  sake  of  completeness,  we  render 
below  specific  expressions  for  a  sphere  composed  of  layered  homogeneous  shells. 

For  convenience,  we  suppoce  that  we  have  always  normalized  nUj  so  that  its  kinetic  energy  is 
equal  to  '/z^u,2.  That  is 


ri 


2  Pt  fKdr  1. 


(78) 


i-1 


Then  the  variation  in  noi,  due  to  perturbations  in  A(,  and  p{  for  /  1 . N  is 

SL a,?)  =  1  [aA  /L dr  f  dpi  /  Mdr  -  /  nw2Kdrl. 

ilL  *,-1  'i-1  'i-1  J 


(79) 


Equation  79  is  an  explicit  representation  of  the  perturbation  in  the  (n,  l)th  squared  eigenfrequency 
as  a  Mnear  combination  of  perturbation  in  the  layered  sphere’s  parameters.  Similarly,  for  Qv  and 
Q'y  ,  we  have 

,1. 


nQ]X  =-y  ^  [<?v'pi  '  2P?  /  Ldr  +  QVsi/'i  / 


(M  -  2L)  dr  .  (80) 
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FORCED  RESPONSE  OF  A  LAYERED  ELASTIC  SPHERE 

We  now  consider  the  response  of ;*  layered  elastic  sphere,  having  a  small  anelasticity,  to  a 
specified  driving  force.  This  species  of  problem  is  classical  (see  Rayleigh  194.))  but  we  give  our 
particular  development  here  because,  as  a  first  approximation  to  our  experiment,  the  specific 
theoretical  results  are  quire  helpful, 

In  the  time  domain,  the  appropriate  field  equation  is 

p  ||  —  L  u  +  L 1  t  I  (81 ) 

where  1.  is  the  elastic  operator  of  eq  II,  I  is  a  spatial  and  temporal  iorce  distribution,  and  L  is  an 
isotropic,  viscous  operator  of  the  form 

L'  =  (A'  +  2/i' )  W  •  -fiVxVx.  (82) 

A1  and  p'  are  taken  to  lie  small  and  constant  (see  eq  68).  Associated  witli  eq  81  are  the  internal 
and  external  boundary  conditions  discussed  earlier. 

We  Laplace  transform  eq  81  and  assume  (for  convenience)  that  both  u  and  fijU  vanish  everywhere 
at  time  t  0.  Then  eq  SI  becomes,  with  p  the  transform  variable. 

(,,p~  -  p  L'  -  Du  f*  (8:i> 

where  u  and  /  refer  now  to  transformed  variables.  We  shall  expand  u  in  the  spheroidal  eigenfunctions 


u 


V 


(84) 


and  the  stun  is  taken  over  all  non-negative  values  of  n  and  /,  and  -  1  <  /  _  I .  The  nu,  are 
orthogonal,  as  a  result  of  the  Hennitiaii  nature  of  L  and  its  boundary  conditions.  The  set  is.  here, 
complete  because  the  set  of  all  free  oscillations  is  complete  (Backus  1967)  and  we  shall  consider 
only  a  class  of  sources  which  cannot  excite  toroidal  modes.  We  may  further  require  the  nuj  to  be 
orthonormal  in  the  sense 


(Pnul-  tiX*  =  /  *V«1 


dv 


~  5nmSlk5js- 


(85) 


If  we  insert  the  expansion  84  into  eq  83,  multiply  on  the  left  by  lnuj. ,  integrate,  and  use  the  re¬ 
sult 


L  u{ 


~P  n1 


we  find,  after  applying  eq  85,  that 
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Since  V  is  an  isotropic  operator,  the  inner  product  (n|uJk.  L'  nujs)  will  vanish  for  k  /  l  or  j  ^  s  as 
a  result  of  the  orthogonality  of  spherical  harmonics.  We  do  not  know  that  this  is  true  for  k  }, 
j  s,  in  /  n  that  is,  we  do  not  know  that  anelaslieity  of  this  sort  does  not  result  in  coupling  be¬ 
tween  overtones  with  the  same  order  1.  For  onr  purposes,  however,  we  are  interested  in  the  effect 
of  the  anelastic  term  when  one  mode  is  strongly  excited  relative  to  other  modes,  i.e.  near  a 
resonance.  This  case  is  approximately  the  same  as  a  response  composed  of  a  single  mode,  and 
we  shall  ignore  possible  overtone  coupling.  Following  eq  70,  we  define  the  scalar 


nDt  (n<  L'  ,/t'- 

Then  we  have 


P  n°l  n‘"l 


(86) 


(87) 


( I  his  does  not  constitute  a  redefinition  of  :|Dj.  The  quantity  defined  above  is  the  same  as  that  given 
by  i  ]  69-70.) 

We  assume,  for  f  (r .  t),  the  form 


f  <r ,  t)  r  S  (r  rs) 


f  sin  (kO  f  ^  0 
1  0  {  •_  0 


(88) 


where  <5 (x)  is  the  three-dimensional  Dirac  delta  function  and  r  is  the  "source"  location.  Then 
eq  87  becomes 


j  *nWyK.<V 

—  - - - - -  (89) 

(p“  *  p  n°\  '  n,‘J^,)<p'',  ’  K') 

n^i  <'„) 's  **’e  val"e  ot  D|e  radial  dis[)laeement  component  at  rs,  the  source  radius.  Since  radial 
displacement  vanishes  identically  for  toroidal  modes,  a  radial  force  source  will  not  excite  toroidal 
oscillations. 

In  order  to  compute  the  time-domain  solution  u(r,  f),  we  have  to  apply  the  inverse  transfoim  to 
eq  89  and  use  the  resultant  „«{(()  in  eq  8d.  The  task  is  straightforward  and  we  mention  here  only  its 
salient,  albeit  predictable,  features.  n«|.  as  given  in  eq  89,  has  four  singularities,  one  at  each  of 
*  is,  one  at  -  14  nD{  +■  i  \Znwf  -  '/*  nDj  and  one  at  -  U  nD{  -  i  \/n<^  -  li  nD|.  The  former  two 
represent  forced  motion  of  the  body  at  the  frequency  #>..  Technically  speaking,  the  forced  component 
of  motion  exists  because  of  the  singularities  in  the  transform  of  the  source  time  function.  Time 
functions  such  as  delta,  ramp  and  step  functions  whose  transforms  lack  singularities  do  not  give 
rise  to  such  forced  motion. 

The  second  set  of  singularities  represents  free  vibration  at  the  body's  various  natural  periods. 
These;  motions  result  from  "turning  on"  the  source  at  t  0.  They  are  damped  exponentially  and 
the  farther  in  time  one  is  from  t  0.  the  weaker  these  motions  become. 

For  our  purposes,  we  require  only  the  forced  motion  and  we  will  neglect  the  remaining  poles. 

The  result  is 
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a  (t) 


nUl{rs]  Y\{0»'  |(„jL>  _  k2)  sin  (Kt)  +  kD  cos  Ut)l 


(90) 


(* 


+  k~D< 


where  the  subscripts  have  been  deleted  from  a,  u>  and  D  for  convenience. 

We  now  suppose  that  the  source  is  located  at  6'  0,  ;  and  that  at  0  =  n,  r 

is  emplaced  a  radial  motion  detector.  For  the  latter  location,  eq  84  becomes 


1  .,</%)(-  1)' 


n  i 


a  +  i 

4  n 


wliere  the  sum  is  taken  only  over  n  and  I  as  a  result  of  the  source  location. 
The  received  power  as  a  function  of  driving  frequency  k  is  given  by 


R2(k)  = 


,/!,(-  lHo/ 


Ov  *> 


r(tu2  - 


2\2  2  2  /n2 

k  )  f  k  u)  W 


i  2  2.2  2  /n2 

(a)  -  K  )  t  K  (O/Q 


(91) 


wliere 


(92) 


nl 


We  have  expressed  R~(k),  the  squared  signal  strength,  in  terms  of  the  source  factors  and  the. 
lumped  quality  factors  nQ,.  If  only  a  single  mode  need  be.  considered,  we  have 


*<„“!>  n^l  ' 

which  relates  response  maxima  to  both  ^  and  nQj. 

It  is  straightforward  to  show  that  when  only  one  mode,  is  significant  and 

R'~(k)  =  i  R2  (_«>,) 

2  n  1 

then  to  first  order  in  nQj 1 ,  «  satisfies 
t  1 

ntJl  o  nwl 


(98) 


(94) 


(95) 


n«l 


Equation  95  is  the  classic  lelation  between  dissipation  and  peak  width,  and  is  useful  only  if  nQj 1 
is  sufficiently  small. 

» 

We  have  supposed  above  that  the  amplitude  of  l  and  the  sensitivity  of  the  radial  displacement 
receiver  were  independent  of  frequency  »<.  If  this  is  not  so  neither  eq  93  nor  eq  95  is  necessarily 
correct,  llowevei.  the  half-power  relation,  eq  95,  will  generally  be  usable  if  source  strength  and 
receiver  sensitivity  do  not  vary  appreciably  over  the  range 
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1,,‘V1  -  f  i mq;‘)| 


g  n^l  '•  n“l 


n**l 


In  antieipat ion  of  experiment,  we  mention  that  if  the  above  influences  are  important  or,  as  is 
more  often  the  ease,  if  the  n«i>[  are  sufficiently  “dense”  on  the  *-axis,  then  neither  eq  92  nor 
eq  95  will  yield  reliable  measures  of  Qjl.  If  this  is  the  ease,  we  will  not  be  able  to  improve  matters 
appreciably  by  making  time  domain  measurements  of  the  decay  of  a  given  mode.  Irregularities  in  the 
“forced"  spectrum  will  appear  in  the  latter  as  distortions  of  a  mode’s  exponential  decay  envelope, 

We  wish  to  make  one  further  point,  regarding  restrictions  on  the  nature  of  the  dissipative 
mechanism,  eq  82.  Suppose  that  the  time  domain  stress-strain  relation  is  expressible  as 

T  =  Te  h  fC(t  -  r)  :  Jre(r)dr,  (96) 

0 

where  T°  >s  the  elastic  stress  defined  in  eq  1,  C  is  a  fourth  order  tensor  function  of  the  time  '‘lag," 
t  -  r,  and  e(r)  is  the  strain  tensor, 

e ( r)  I  (Vu  .  UV).  (9") 

-  2 

We  suppose  both  stress  and  strain  to  vanish  for  negative  times.  Laplace  transforming  eq  96  and 
supposing  0(0)  to  vanish  gives 

T  T°  ,  p  C  (p)  :  e(p).  (98) 

When  C  is  an  isotropic  tensor,  the  use  of  the  stress  defined  in  eq  9S  ui  the  transformed  equations 
of  motion  (eq  2)  will  yield  exactly  operators  of  the  form  L  and  L1  of  eq  81  and  82.  These  differ  from 
our  usage  only  in  that  A1  and  p'  may  now  depend  functionally  on  p,  the  transform  variable.  The 
quantity  t)D|  defined  by  eq  86  or  eq  69-70  would,  in  turn,  depend  upon  p.  If  ,,0,  (or  A'  and  /t’)  does 
not  vary  appreciably  over  a  frequency  span  of  a  “few  peak-widths,"  the  results  derived  above  should 
be  usefully  accurate.  (We  must  still  have  A'  and  /«',  or  nD,,  small  since  we  have  used  a  perturbation 
method  to  include  attenuation.) 


THE  INVERSE  PROBLEM 

In  two  of  the  preceding  sections  we  developed  the  "forward''  problem  for  a  layered  elastic 
sphere.  The  forward  problem  consists  of  the  generation  of  the  eigenvalue  spectrum  associated 
with  a  given  model.  In  this  section  we  consider  the  inverse  problem  of  utilizing  a  measured 
eigenvalue  spectrum  (which  will,  inevitably,  be  incomplete)  to  infer  the  properties  of  the  materials 
of  which  the  sphere  is  composed. 

Let  M  be  the  space  of  all  layered  elastic  spheres  having  N  layers  delimited  by  the  points 
lr0,  r1(  . . rN  t  where,  as  before,  rQ  =  0  and  rN  is  the  sphere’s  radius.  Then  all  models  in  Af  have 
a  common  geometry  but  differ  in  the  elastic  properties  (including  density)  of  their  component  shells. 
M,  then,  is  a  space  of  dimension  3 N  (since  each  shell  has  three  distinct  properties)  and  we  may 
represent  a  given  model  by  m,  a  vector  of  dimension  3 N.  We  further  limit  M  by  requiring  that  it 
encompass  only  physically  realizable  models.  A  model  is  said  to  be  physically  realizable  if  each 
shell’s  properties  satisfy 
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p  s  0. 
P  1  0. 


and 


•) 


(99a) 

(991)) 


(99c) 


whore  both  equalities  99b  and  99c  arc  not  simultaneously  true.  The  latter  two  constraints  merely 
express  the  condition  that  an  elastic  material  be  thermodynamically  stable  (Font;  1 965). 

Let  n<a|(m)  be  the  n,fl  overtone  of  the  spheroidal  mode  of  degree  /  associated  with  the  model  m. 
( Both  n  and  l  range  over  ttie  non-negative  integers.)  We  cannot  express  u».  (m)  in  closed  analytic 
form  but  we  can.  through  techniques  previously  discussed,  generate  it  numerically.  We  can  now 
formulate  ttie  inverse  problem  in  the  following  manner  (Backus  and  Colbert  1967). 

Let  w!1.,  i  1 . K  he  observed  resonant  frequencies  associated  with  particular  modes  of 

oscillation.  We  wish  to  determine  a  model  m,  satisfying 


in  It 


o 

m“’li 


1 . K. 


(100) 


As  Backus  and  Gilbert  (1967)  have  pointed  out,  we  do  not  know,  a  priori,  if  the  set  of  solutions 
to  eq  100  is  empty,  has  a  single  member,  or  is  a  subspace  of  M  of  one  or  more  dimensions. 

We  do  not  know  a  direct  procedure  for  solving  eq  100  lor  one  or  more  models  m.  We  resort  here 
to  iterative  methods  for,  hopefully,  generating  successively  improved  approximations  for  m.  For 
convenience  we  rewrite  eq  100  as 


D((m)  -  D“  j  1 . K  (101) 

where  the  are  data  and  the  0,(m)  are  data  functions,  D,( m)  is  a  scalar- valued  function  whose 
domain  is  the  3N- dimensional  vector  space  M  and  whose  value  is  the  angular  frequency  of  free 
oscillation  of  the  n{ft  overtone  of  degree  I  .  Let  m0  be  some  model  which  we  believe  to  lie  near  m. 
We  wish  to  find  some  perturbation,  5m,  m  mn,  such  that  m°  *  5m  more  nearly  satisfies  eq  101. 
(m°  and  m°  *  5m  must  both  lie  in  M  but  5m  alone  need  not.)  We  wish  to  have 


D{  (m°  *  5m)  =  Of  i  1,  . . .  K. 

Expanding  D^( m)  m  a  Taylor  series  gives 

3  N  \0D 

DAm°  *  5m)  =  DAm°)  •  V*  — 

1  '  i—i  15m, 


(102) 


j-1  *  J 

Then,  to  first  order  in  5m j,  we  wish  5m  to  satisfy 


8m,  *  OCSm  ") 
o  J 


i  -  1. 


(103) 


h  K 


5m,  0['  /^(m0)  i  l,  ...  K. 


(101) 


1 
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II  |Am|  is  sufficiently  small  we  may  ex|ieot  that  m‘  m11  •  Am  will  mufc  nearly  satisfy  eq  tot 
t liun  m°  d:J.  As  a  measure  of  a  model's  suitability.  we  may  define 


F.quations  101  constitute  a  K  *  :i N  set  of  linear  equations  in  the  components  ot  Om  We  may 
not,  in  general,  expect  to  find  Am  exactly  satisfying  eq  101  tor  all  |mssible  cases.  If  the  rank  of  the 
system  does  not  exceed  UN.  such  a  Am  exists  tint  is  not  necessarily  unique.  If  the  rank  exceeds 
3 N,  it  does  no,'  exist. 

Vario  is  methods  of  solution  have  lieen  applied  to  the  system  101  (Backus  and  flilhert  1907, 
Anderson  and  Smith  1908.  Smith  and  Franklin  1909.  Jordan  and  Franklin,  manuscript,  1971).  We 
adopt  here  a  general  technique  promised  hv  Franklin  (inqnddished  inaiiuscript.  1969). 

We  rewrite  eq  104  mote  compactly  as 


/l(m°l  •  Am  Rim11) 


(100) 


where  A(m°)  is  the  matrix  whose  elements  a  are 


atj(") 


AD, 
Am , 


m 


(107) 


and  R(m°)  is  the  vector  of  data  residuals  whose  i'h  component  is  Du  -  D^m0).  We  now  regard 
eq  106  as  a  linear  relation  between  tlnee  stochastic  processes  a  signal  pioeess,  a  data  process 
and  a  noise  process.  Am  is  a  sample  of  the  signal  process,  and  R(m°)  is  the  sun  of  a  sample  of 
the  data  process  plus  a  sample  of  the  noise  process.  Kadi  process  is  taken  to  have  zero 
expectation. 

The  use  of  stochastic  techniques  to  solve  eq  106  is  based,  in  part,  upon  contemplation  of  some 
of  the  potential  sources  of  error  entering  into  the  relation.  The  measured  spectral  values  D®  are 
contaminated  by  measurement  error  and  possible  mode  misidcntificurioii.  The  physical  system  upon 
which  measurements  are  made  may  deviate  from  t Ire  class  of  models  M  in  which  m  must  lie.  It  is 
possible  that  there  is  no  model  m  in  M  that  would  then  satisfy  eq  101. 

Let  Rm  denote  the  autocorrelation  associated  with  the  signal  process  Am  ,  and  R°  denote  the 
autocorrelation  operator  associated  with  the  noise  process.  Km  is  a  3N  v  3JV  square  matrix  whose 
(i,  j)th  component  is  the  expectation  of  the  product  A/rySnij.  jP  is  defined  analogously.  The  best 
linear  estimate  for  Am  is  given  by  (Franklin  1969): 


Am  -  Rm  ■  /lT(m°)  •  U<m°)  •  Rm  •  4T(m°>  ♦  K0!"1  •  R(»°).  (108) 

A  solution,  Am,  can  be  guaranteed  to  exist  if  R(l  is  a  positive  definite  matrix. 

The  above  method  for  computing  Am  was  chosen,  in  preference  to  such  techniques  as  least- 
squares  solutions  because  experience  has  indicated  that  the  solution  108  is  typically  smaller  and 
more  stable  than  that  provided  by  better  known  methods.  Since  eq  108  is  used  only  to  compute  the 
perturbation  Am,  which  we  use  to  iteratively  improve  m(\  it  does  not  follow  that  the  enors  in  our 
"final"  fit  will  correspond  to  any  of  the  components  of  K(l  or  that  the  “distance"  between  our  initial 
and  final  models  will  he  strongly  related  to  Rm. 
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In  practice,  we  take  if  to  be  a  K  •  K  diagonal  matrix  whose  i"'  diagonal  element  is  the  square 
ot  our  estimate  of  the  uncertainly  of  the  i'h  datum.  R"'  is  also  taken  to  he  a  diagonal  niatiix  whose 
i,h  flement  is  taken  equal  to  •"(if.  where  t  is  a  small  number  (we  used  0.2  to  0.01),  and  </,  is  the 
component  of  m.  Thus,  we  “suggest"  to  eq  108  that  changes  in  a  given  component  or  m  should 
be  of  the  older  of  a  fraction  of  its  value.  (Life  is  not  quite  that  simple  since,  as  Franklin  (1069) 
points  out.  it  ts  only  the  ratio  of  signal  to  noise  that  counts,  but  the  values  given  above  tend  to 
work  well  m  practice.) 

It  is  also  helpful  to  have  some  measure  of  the  uncertainty  associated  with  results  so  obtained. 
Such  estimates  must  rclv.  of  course,  on  the  statistical  properties  we  assign  to  errors  m  the  data 
and  errors  in  the  model  and  Hie  fashion  in  which  one  relates  to  the  other.  Forgoing  this,  we  adopt 
here  an  estimate  which  is  simple,  and  perhaps  crude,  in  the  extreme. 

We  define  the  r.m.s.  absolute  error  in  oigenlrequency  i,  as 

L  i  K>).  <109> 

'  E  ,  ,  1 

Associated  with  each  component  hi,  of  m  is  an  r.m.s.  sensitivity  defined  by 

,r  L  1  aj>).  <1,0> 

We  take  as  a  measure  of  the  uncertainty  m  If  the  data  depend  only  weakly  upon  m,.  then 
„  will  lie  small  and  .f  </,  large.  Similarly,  tf  the  data  depend  strongly  on  rf,  will  he.  (relatively) 
large.  We  cannot  offer  a  more  quantitative  justification  for  this  method. 

EXPERIMENTAL  AND  NUMERICAL  RESULTS 


Illustrative  applications 

hi  order  to  render  more  specific  the  developments  of  preceding  sections,  we  give,  first,  the  re¬ 
sults  of  two  experiments  and  various  aspects  of  their  interpretation.  The  first  sample  discussed 
here  was  a  l-in.-diutu  sphere  of  Lueite.  formed  by  machining  from  a  section  of  bar  stock. 

Figure  1  shows  the  observed  forced  spectrum  of  this  sample  as  recorded  by  a  swept -treqiieucy 
analyzer.  The  sample  was  moderately  clamped  between  a  quartz  transducer  and  a  PZT-I  transducer , 
which  happened  to  be  available,  and  was  in  a  refrigerated  cabinet  at  about  10  C.  ’lhe  numbers  are 
peak  frequencies  in  kilohertz,  as  measured  with  a  counter,  and  the  mode,  assignments  were  made  on 
the  basis  of  considerations  ontliued  below. 

From  past  exjrerience .  it  was  possible  to  readily  identify  the  “main  sequence"  of  modes  nS.,. 

S . all  of  which  have  substantial  maxima  aud  which,  for  1  in  excess  of  2.  tend  to  be  more  oi  less 

evenly  spaced.  Since,  as  we  shall  show,  the  projmrties  of  this  sequence  arc  generally  dominated  by 
shear  velocity  over  regions  of  Poisson’s  ratio  of  common  interest,  it  was  desirable  io  identity  at 
least  one  of  0S„  or  0S,.  both  of  which  are  significantly  Influenced  by  compression;.!  velocity,  hi 
identifying  0S0,  we  applied  two  useful  guides,  gleaned  from  contemplation  of  eq  91. 

1.  Adjacent  modes  both  having  /  even  or  I  odd  will  interfere  destructively. 

2.  Adjacent  modes  not  both  having  the  same  parity  will  interfere  constructively. 

These  guides  applv  to  frequencies  lying  between  the  two  appropriate  eigeufrequencies.  In  light  of 
these,  we  interpret  the  sharp  asymmetry  of  (1S„  as  being  a  resell  of  destructive  interference  with 

„S4  (Mi  the  left. 
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Figure  1.  Observed  forced  response  spectrum  of  a  4-in.-diam  Lucite  sphere  at  - 40°C .  Ordinate  is 
linear  and  arbitrary.  Mode  assignments  and  peak  frequencies  are  shown  in  kilohertz. 


From  these  data,  we  estimated  that  Vp  was  about  2.84  km/sec  and  l's  was  about  1.44  km/sec. 
Using  these  values  we  carried  out  forward  calculations  for  all  spheroidal  modes  having  frequencies 
less  than  105  Hz,  the  results  of  which  are  shown  in  Figures  2  and  ft.  Appendix  C  discusses  some 
numerical  aspects  of  the  forward  calculation. 

Figure  2  depicts  the  source  factors,  nA1,  for  the  fundamental  and  first  four  overtones.  These 
results  confirm  the  identity  of  the  main  sequence  and  qSq,  and  permit  the  identification  of  a  number 
of  additional  modes.  Most  of  these  are  fairly  deeply  buried  in  the  "noise”  and  are  useless  as  data 
but  they  do  help  confirm  our  identifications.  Prior  to  the  use  of  these  source  factors,  and  the  inter¬ 
ference  rules,  we  often  found  that  our  results  were  marred  by  uncertainty  about  mode  identification. 
The  above  aids  have  greatly  improved  this  aspect  of  the  technique. 

Figure  3  shows  the  quantity 


K  + 

where  Qa  and  Q ^  are  the  weights  computed  from  eq  77  to  average  the  invers^Q’s  of  compressional 
and  shear  waves  to  arrive  at  a  mode's  lumped  inverse  Q.  (Actually,  Qa  t  Qp  1,  but  t lie  above 
form  is  more  indicative.)  In  addition  to  being  inherently  significant,  this  quantity  is  a  convenient, 
normalized  index  of  the  partition  of  energy  into  compressional  and  shear  waves  and  thus  of  the 
properties  “controlling”  a  given  mode.  (The  two  families  nS0  and  nS,  have  been  joined  by  separate 
curves.)  Figure  3  suggests  that  the  main  sequence  is  dominated  by  shear  velocity  and  almost  uni¬ 
formly  so.  We  should  not  expect,  then,  to  get  good  values  for  compressional  velocity  from  this  set 
alone.  The  set  nSj  for  n  2,5,8  etc.  is  controlled  by  V  but  does  not  usually  show  up  well  on 
observed  spectra.  Experience  indicates  that  0S0,  and  to  a  lesser  extent  qSj,  provide  the  most 
commonly  observed  control  on  V  .  Other  modes  are  useful  only  when  Q  is  sufficiently  high  that 
they  appear  clearly. 


Figure  2.  Source  factors,  j4j,  for  spheroidal  modes  of  a  4-iu. -diuni  Lucite  sphere  at  40  ‘C. 
Ordinate  is  logarithmic.  Arrows  for  I  0  and  I  1  indicate  the  overtone  assignments  for 
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Figure  4.  Portion  of  the  theoretical  forced  response  spectrum  of  a  4-in.-diam  Lucite 
sphere  assuming  Qa  =  5 00.  The  peak  is  0SQ  and  is  asymmetric  due  to  de¬ 

structive  interference  on  the  left  with  0$4. 

In  this  instance,  as  Figure  1  shows,  Q  is  so  low  that  we  cannot  set  a  baseline  in  order  to  esti¬ 
mate  peak  widths.  For  higher  Q  materials,  we  are  sometimes  able  to  measure  mai  l  sequence  attenua¬ 
tions  from  peak  widths.  We  are  seldom  able  to  get  reliable  (?’ s  for  qSq  or  other  compression- 
dominated  modes.  Similar  difficulties  were  reported  by  Fraser  (1970)  in  vitreous  silica,  and  it  does 
not  seem  likely  that  we  will  be  able  to  get  reliable  compress ional  attenuation  data  in  any  simple 
fashion, 

This  gloomy  statement  is  bolstered  by  theoretical  computations  of  the  forced  response  using 
eq  91  and  assumed  values  of  Q(l  and  Q^.  Figure  ■!  shows  the  theoretical  j)S4-()S0  interaction  when 
Qa  and  Q ^  are  both  equal  to  500.  (These  results  are  based  on  calculations  for  the  -l-in.  Lucite 
sphere.)  In  addition  to  the  resemblance  of  Figure  1  to  Figure  1,  we  note  that  „S0  still  possesses 
substantial  asymmetry  even  though  the  two  peaks  are  (in  this  case)  about  32  peak-widths  apart. 

Figure  5  depicts  the  interaction  of  ,S4,  L,Sj  and  2S 2  with  0S?  on  the  left  and  ().Sg  on  the  right. 
Calculations  for  two  values  of  Q(l  ~  Q ^  are  shown.  The  asymmetry  of  each  of  5S4  and  „S(  is 
clearly  due  to  interaction  with  the  more  distant  main  sequence  inodes  rather  than  with  each  other. 

The  appearance  of  all  three  modes  is  dominated  bv  their  positions  on  the  "tail  of  0S7.  The  figure 
also  depicts  the  rapid  degradation  of  low  amplitude  modes  with  decreasing  Q.  The  maxima  associated 
with  ,S4  and  2S2  become  increasingly  blurred  as  Q  decreases.  In  the  case  of  Lucite,  which  has  a 
compressional  Q  of  about  70,  Figure  1  shows  that  only  „S(  jiersists  as  a  noticeable  maximum. 

We  attempted,  next,  to  "invert"  the  data  obtained  from  the  main  sequence,  plus  ()S0  in  three 
different  ways.  The  results  of  this  endeavor  are  given  in  Table  1. 

The  first  result,  Inverse  1,  was  obtained  by  using  all  the  data  in  the  automated  algorithm.  The 
"fit"  has  an  r.m.s.  relative  error  of  l"i,  due  largely  to  the  2.6"',  error  in  matching  0S„.  Inverse  It 
was  obtained  by  manually  matching  ()S0  and  This  result  has  the  same  r.m.s,  refative  error 
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Figure  5.  Portion  of  the  theoretical  forced  response  spectrum  of  a  4-in.-dinm  Lucite  sphere  for 
two  assumed  values  of  Qu  Qq.  The  figure  depicts  the  effects  of  interference  on  0S4,  2Sj, 

and  gS>)> 


Table  I.  Comparison  of  velocity  estimates. 


Mode  Observed  Inverse  I  Inverse  2  Inverse  3 


kHz 

kHz 

kHz 

( 

kHz 

0^0 

24.27 

24.41 

-r>.8 

V 

10'1 

24.29 

-8.2  > 

10"1 

24.27 

1.28 

* 

10' 

<A 

11.54 

11.84 

-2.6 

■% 

10'1 

1 1.92 

-0.3  x 

10“ 

11.91 

-3.2 

y 

10' 

oS, 

17.78 

17.64 

7.9 

V 

10'’ 

17.77 

5.6  v 

10'* 

17.75 

1.65 

10' 

22.78 

22.63 

6.6 

10' 1 

!  2.79 

—4.3  • 

10“ 

22. 78 

3.99 

X 

10' 

<A 

27.43 

27.30 

4.7 

< 

10'' 

27.50 

-2.5  \ 

10'' 

27.47 

-1  44 

« 

10' 

qS* 

32.01 

31.82 

5.9 

y 

10' 1 

32.04 

-0.4  v 

10“ 

32.01 

-41.0 

10' 

A 

36.13 

36.25 

4.9 

y 

10' ' 

36.50 

-1.9  s 

10“ 

36.46 

-8.7 

• 

10' 

oS, 

40.87 

40.62 

6. 1 

- 

10* 

40.90 

-7.3  v 

10“ 

40.86 

3.2 

V 

10' 

oS* 

45.26 

44.95 

6.8 

*. 

10' 1 

45.26 

0 

45.21 

1.0 

X 

10' 

0^10 

41).  51 

49.26 

r».o 

10' 1 

49.60 

-1.8  v 

10“ 

49.54 

-6.9 

y 

10' 

t 

Obs-Comp 

KMSRK 

1.0 

. 

1<>“ 

KMSRK 

1.0  y 

10'J 

RMSKK 

1.0 

V 

10' 

Otis 


Observed  5.08-cm-radius  Lucite  sphere  at  -40JC 

Inverse  1  Cotrputed  using  all  data  shown 
V  284262 

VH  142577 


Inverse  2  Kstimated  for  fit  to  0S0  and  0S, 
yp  =  285864 
K,  143650 

Inverse  3  Computed  using  all  but  0S2 
V  =  283566 
V's  =  143504 
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despite  an  increase  in  the  error  for  0S.„  and  further  suggests  that  0SL>  is  behaving  anomalously.  The 
third  result.  Inverse  111.  was  obtained  by  using  all  the  data  except  0S„.  and  still  has  an  r.m.s.  relative 
error  of  1%.  The  constancy  of  this  number  is.  so  far  as  we  know,  fortuitous.  However,  if  we  omit 
S.,  entirely,  the  error  drops  to  0.09"u.  It  is  not.  in  our  experience,  unusual  for  the  mode  0S.,  to  he 
inconsistent  witli  remaining  modes.  We  cannot  offer  an  explanation  of  this  observation  at  this  time 
but  do  recommend  that,  when  sufficient  alternative  data  are  available,  it  not  be  utilized.  We  may 
observe,  however,  that  of  the  suite  of  surface  wave  modes.  ()S.,  etc..  ,,S„  does  possess  the  small¬ 
est  surface  wave  number.  (I  i  H)/a. 

Table  11  shows  the  result  of  an  inversion  ol  data  from  a  l-in.  CJK-12f)  fused  quartz  sphere  at 
19‘  G.  The  error  is  0.0.T»  and  „.S„  does  not  show  anomalous  behavior.  These  results  arc  unusually 
Hood.  We  note  t hat  the  relative  errors  surest  that  computed  tioqueitcics  tor  the  main  sequence  are 
deviating  systematically  above  observed  tioqiienrics.  ’once  modes  ol  the  form  (,St  are  known,  lor 
increasing  I,  to  become  nieieusttigly  concentrated  neat  the  sail  ace.  wo  may  speculate  that  such 
systematic  errors  imply  a  lower-velocity  snrttclal  icgintt.  Other  candidates  are  gcomcti ic.il 
irregularities.  innueiiee  of  the  transdueeis.  etc.  WY  tnti.shi.  <-  this  point  because  Fraser  (1970)  has 
shown  that  vacuum  heating  of  some  vitieous  silicas  sttongly  Intluenees  the  attenuation  ol  various 
torsional  (shear)  modes. 


Table  II.  Inversion  of  data  from  a  ♦- In. -d  fa  meter  sample  of  GK-125  fused  quartz  at  19°C. 
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Interferometric  techniques  such  as  this  arc  particularly  useful  for  gauging  tl'c  influence  of 
variables  such  as  temperature  on  velocities.  It  is  often  possible  to  "track  variations  in  velocities 
which  are  smaller  than  the  errors  of  measurement.  Kacli  additional  set  of  data  is  easily  interpreted 
as  a  perturbation  of  the  previous  set  and  does  not  require  additional  mode  identification.  1  inure  6 
shows  results  obtained  from  Limit e  and  fused  quartz  at  three  different  temperatures.  The  data  are 


shown  as  variations  relative  to  values  measured  at  21  C, 
in  the  cabinet  at  *  1°C  and  these  are  depicted  as  horizon! 


We  lave  estimated  temperature  variations 
al  error  bars,  The  linearity  of  the  tempera¬ 


ture  dependence  of  elastic  velocities  in  Lueitc  is  striking.  The  slight  negative  temperature 
dependence  of  fused  quartz  is  consistent  with  results  described  by  Mason  (19.r>8).  (These  data  arc 
shown  purely  for  purposes  of  demonstration.) 
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Figure  6.  Inferred  relative  variations  of  compressional 
and  shear  velocities  in  4-in.-diam  spheres  of  Lucile  and 
GE-125  fused  quartz,  with  temperature.  Results  are 
shown  as  a  relative  change,  for  each  quantity,  from  its 
value  at  21  °C, 


Some  geawal  results  aid  a  graphical  inversioa  method 

Figure  El  (Appendix  E)  depicts  the  variation  of  the  dimensionless  frequency  or  various 
spheroidal  modes  of  a  homogeneous  sphere  as  a  function  of  the  ratio  of  shear  and  compressional 
velocities.  The  dimensionless  frequency  is  related  to  the  actual  frequency  „f,  (we  have  used 

Hertz  here)  by 


n't 


In  addition.  Poissons  ratio  is  shown  for  those  values  of  V/V p  for  which  it  is  positive. 

The  loci  of  all  modes  shown,  except  0S0,  appear  to  converge  on  thefWipina^ 

-  -  ‘  remains  substantially  linear  for  Vs/Pp  as  great  as  0.625.  Clearly  modes 


functions  Ol  KVi'  .  remains  ™  s  p  w 

which  behave  in  tfis  fashion  are  virtually  independent  of  compressional  velocity  and  Figure  El  con¬ 
firms  our  earlier  statements  about  the  sensitivity  of  various  modes  to  compressional  and  shear 
velocities.  The  left-hand  side  of  the  figure  would  have  been  more  manageable  if  1%  had  been  used 
to  nondimensionalize  frequency  instead  of  V  We  took  the  point  of  view,  however,  that  in  experi¬ 
mental  work  one  is  more  likely,  if  any  prior  data  are  available,  to  have  an  estimate  for  Vp. 
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On  the  other  extreme,  calculations  were  carried  to  the  "Limit  of  Elastic  Stability,”  indicated 
in  the  figure.  This  limit  is  imposed  by  the  condition  of ;;  non-negative  compressibility  and 
corresponds  to  a  Poisson's  ratio  of  -1.  So  far  as  we  are  aware,  any  requirement  of  a  non-negative 
Poisson’s  ratio  is  based  strictly  on  intuitive  and  empirical  grounds. 

Figure  E2  is  a  plot  of  the  ratio  of  frequency  for  several  pairs  of  modes  (dimensionless  or  other¬ 
wise)  as  a  function  of  ^s/Vp.  The  ratios  oW:p  ai|(*  :lrt!  SI'0W11  using  the  scale 

on  the  left,  and  qL/q!<>  is  shown  using  the  greatly  expanded  scab*  on  the  right.  Figure  ED  implies 
that  the  pair  0f2  and  0/2  will  not  serve  well  to  determine  l'B/P  whereas  „/„  and  ()f;i  or  ()/2  will. 
The  pair  0/ 1  and  QL,  is  of  ques  lonahle  value. 

The  plots  are  useful  ns  a  means  of  graphical  inversion,  as  we  shall  now  show  from  the  data  of 
Table  Ill.  For  the  4-in.  fused  quartz  specimen  we  have 


and 


0*0 

44555 

0^2 

30898 

0^0 

44555 

0^3 

45580 

(/.I 

45580 

1 .4120 


0.977!) 


0/2  30898 


=  1.47517 


which  imply  (from  Fig.  E2) 

V's 

—  =  0.651,  0.633  and  0.631. 
I> 


For  V  /V  -  0.632,  Figure  El  gives 

x  r 


0^0  rN 


0.3807 


or 


591535 


and 


Vs  =  375746 

which  differ  about  0  2%  from  the  values  determined  by  automated  inversion. 

This  method  should  provide  usefully  accurate  results  when  individual  data  are  of  a  fairly  high 
quality  -  as  is  often  the  case  for  isotropic,  high-Q  materials.  When  only  data  of  lower  quality  are 
available,  or  greater  precision  would  be  useful,  we  must  use  some  sort  of  automated  algorithm. 
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Table  DI. 


Forward  computation  for  1’  593980  cm  sec 

37  5540  cm  sec 
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Automated  Inversion 

Tlie  inversion  of  spectral  data  to  yield  elastic  velocities  is,  e\cit  in  onr  relatively  simple  case, 
a  fundamentally  nonlinear  process.  Consequently  it  is  difficult  or  impossible  to  provide  the  proofs 
of  uniqueness  commonly  available  for  linear  problems.  The  more  general  problem  of  uniqueness  when 
elastic  properties  are  functions  of  radius  has  boon  extensively  discussed  (see.  for  example,  Backus 
and  Gilbert  1967)  but  the  results  are  not  of  much  use  to  ns  here. 

The  monotonio  dependence  of  the  ratio  of  0f0  to  0L„  and  of  the  dimensionless  frequency  0(q,  on 
V'/V  suggests  that,  within  the  range  of  shown  in  Figure  E2,  inversion  is  unique  if  the  proper 

data  are  used. 

In  a  realistic  application,  however,  wc  must  account  for  several  additional  factors.  The  follow¬ 
ing  effects  may  be  important 

1.  Measurement  errors  or  contamination  of  the  data  by  asphericitv,  anisotropy,  inhomogeneity, 
a  finite  Q,  and  external  influences. 

2.  Mode  misidentification. 

3.  Numerical  instability  in  the  computing  processes.  The  result  of  either  1.  or  2.  may  be  to 
produce  a  set  of  data  which  no  homogeneous  sphere  will  "fit,'’  The  influence  of  such  errors  will 

be  heightened  if  we  attempt  to  use  contaminated  data  to  determine  a  parameter  which  is  only  a  slowly 
varying  function  of  these  data.  Such  a  problem  is  said  to  be  "ill-posed  ’  and  the  effects  of  small 
errors  in  "ill-posed"  computations  arc  often  quite  large  (see,  as  an  example,  Franklin  1968). 

While  we  have  not  devoted  extensive  attention  to  this  problem,  we  have  performed  a  few 
exploratory  calculations.  Table  111  shows  the  result  of  two  inversions  using  artificial  data.  The 
starting  models  were  chosen  to  be  about  20%  "away"  from  the  known  correct  model.  "E"  denotes 
the  r.m.s.  relative  error  of  the  data  at  the  end  of  each  pass.  The  process  is  rapidly  convergent.  We 
have  also  performed  similar  calculations  with  data  to  which  a  1%  random  error  was  added.  Such 
inversions  no  longer  gave  the  correct  answer  but,  in  the  few  cases  we  tried,  the  answer  remained 
independent  of  the  starting  model. 
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It  is  easy  to  conceive  of  various  Monte  Carlo  experiments,  requiring  the  usual  spectacular 
amounts  of  computing  time,  which  would  help  delimit,  for  a  given  data  set,  the  error  level  at  which 
inversion  becomes  “incoherent."  Because  of  the  large  multiplicity  of  parameters  required  to  cover 
just  typical  experimental  cases,  we  do  not  presently  think  it  worthwhile  to  do  this.  As  a  matter  of 
experience,  we  have  never,  so  far  as  we  can  tell,  encountered  any  pathological  results. 


Errors 

We  will  briefly  discuss  only  one  potential  source  of  error;  the  effect  of  mechanical  loading  by 
transducers  or  support  devices.  Appendix  D  describes  a  computation  which  purports  to  estimate 
the  perturbation  in  eigenfrequency  induced  by  clamping,  with  a  given  force,  a  sphere  against  a 
halfspace  of  known  elastic  properties. 

The  errors  so  deduced  do  not  exceed  a  few  parts  in  104  for  a  fairly  realistic  range  of  properties. 
We  believe  these  results  lend  support  to  claims  by  Fraser  and  LeCraw  (1964)  and  Birch  (personal 
communication)  that  this  method  is  relatively  free  from  errors  arising  from  such  sources. 
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APPENDIX  A.  ANALYTIC  SOLUTIONS  OF  THE  TRANSFORMED  WAVE  EQUATION 

We  represent  u  by 

u  =  yi/»  +  V  x  ~ffX  +  V  x  V  xTro  (Al) 

where  i/»,  x  and  a  are  scalar  fields.  The  representation  A1  was  chosen  because  it  is  "natural”  to 
the  field  equations  (eq  3). 

It  is  helpful  to  develop  the  following  useful  expressions 

V  •  u  =  v2<A. 

and 

yxyxu  =  yxVxyx  ("X)  +  V  x  V  x  V  x  Vx  (Fra)  = 

=  V  x  lv(V  •  Try)  -  V2(rrx)l  + 

+  v  x  V  x  lv(  V  •  Tra)  -  V2(rra)l . 

We  expand  V2  (f 1 X )  as 

=  Tr  V2X  +  Vl2y! 

as  may  easily  be  shown  by  expanding.  So, 

y  x  vxu=-vxFrv2X:-Vxvxr"r  V2a. 

If  we  insert  these  into  the  field  equation  (eq  3)  and  regroup  terms,  we  find 

VI  pd^ifj  -  (A  +  2/i)V2  <AI  +  V  x  Tt  ipdfx  -  ^V2xl  + 

+  y  x  y  x  7r  I  p  <?2  a  -  p  y2  a  i  =  0. 

In  order  to  ensure  that  eq  3  is  satisfied,  it  is  sufficient  that  i/< ,  x  and  a  be  solutions  to 

pd *if,  -  (A  +  2 p)v2<A  =  0, 


(A2a) 
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P#2X  ~  PV2X  =  0,  (A2b) 

and 

pd*o  -  p v2a  =  0.  (A2c) 

We  have  not  shown  that  all  solutions  to  eq  3  can  be  expressed  in  terms  of  functions  satisfying 
etj  A2  and  Al.  These  solutions  are,  however,  known  to  be  complete  (Sternberg  1960). 

We  now  Fourier  transform  the  system  and  introduce  the  expansions 

ip  (r,  6,  tf>,  u)  =  v  v  (r>  w)  ym  (0,  0),  (A3a) 

1=0  m=-l 

xr<f.  e,  <f>.  u>)  =  k  x?('.u)YT  <0.0).  (A3b) 

1=1  in  =— 1 
and 

or(r,  0,  <f>,  <j)  =  £  1*  o™(r,  w)  V'|n  (0,  0).  (A3c) 

1=1  m=-l 

The  terms  of  degree  1  =  0  in  the  expansions  for  \  and  °  have  been  omitted  since  they  do  not 
contribute  to  the  displacement  field. 

The  expansions  A3a-A3c  are  inserted  into  the  transformed  versions  of  A2a-A2c.  We  make  use 
of  eq  6  and  14  to  simplify  the  result.  If  (  is  some  scalar  field,  then  by  eq  6 

Vf  =  70/  +  Vj  (r-1 0 
and 


V2f  =  d2(  +  j  0/  +  7*^*  1 


.-2  —  2 


(A  4) 


The  resulting  expressions  are  multiplied  by  ( 0 ,  <£)  sin  0  and  integrated  over  0  and  .  We  appeal 

to  the  orthogonality  relation  12.  We  then  have 


l  f  r  r  A  +  2,  r2  f 


(r.  w)  =  0, 


>)  =  0, 


(A5a) 


(A  5b) 


and 
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1(1 ,  lij.j.ip  it,  „>)  -  o. 


Each  of  the  operators  in  brackets  is  some  form  of  the  spherical  Bessel's  operator, 
are  well-known  and  they  are 


0™  =  A™  Uo)  ix(kr)  t  B"1  (.*>)  y,  (kr). 
X™  =  C“(o>)),lyr)  4  DfUy,  (yf). 


(A5e) 


Its  solutions 


(A6a) 

(A6b) 


and 


a*  =  U)/,(yr)  4  Ff  («.)y,(yr). 


(ABc) 


where 


(A7) 


(A8) 


We  Ml  now  to  relate  tt.V  and  »  to  *  y  and  ...  To  do  thin,  we  will , narrate  .<1  Al  to 
resemble  eq  4.  For  convenience  we  will  drop  subscripts  and  superscripts. 

We  note  immediately  that 
V0  =  7*3r0  4- 


Also, 

V  x  Tr  x  =  r  XV  x  r-rx  V(rX) 

=  -T  X  VO"*) 

=  -  r~ X  V!X. 

The  third  term  can  be  expanded  as 

y  x  v  xTrt;  =  V(V  •  rro)  -V’frra) 

=  -Trv2"  -  V(2a)  +  v  lr-s «9r  (r3rr 
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V  *  v  rr«»  r*l  -  r 

*  lr":,^r(r:,n)  -  2r~li»  I 

r*l  r~^'n  I  •  +  <>t<i  I 


.(-HI  #  1)  \  _  v| 

| - ,t  v  .  ^jlr  (r  fi ) ! . 


Collecting  these  we  lutve 


u  IU ■  |  >  VjirV*  *  r_l  <if  <rw  1 1  -  r  V,.y. 


Therefore 


t/jn  ,  6rtf‘  -  «T, 


V™  n  r->  [ C,®  .  ,lf  (rr/f  )|. 


and 


u/ni  .  ,m 

wi  \t 


(A9;t) 


(A91>) 


(A9c) 


APPENDIX  B.  THE  TRACTION  ON  A  SPHERICAL  SURFACE 

Let  T  be  the  elastic  stress  tensor  resulting  from  a  displacement  field  u.  T  is  given  by 

T  =  A(V  •  u) I  +  f«(V u  +  uv)  (B1) 

where  uV  is  simply  the  transpose  of  Vu.  The  force  (not  stress)  acting  across  a  surface  whose 
unit  normal  islTis  given  by  7f  •  T ,  which  is  a  vector  quantity.  We  represent  the  force  acting 
across  a  surface  whose  normal  is  the  radius  vectorTby 

r-  T  =  TP  +  Vj Q  -  r* x  V^.  (B2> 

Our  problem  is  to  relate  the  scalars  P,  Q  and  R  to  the  scalars  U,  V  and  W  which  characterize  the 
displacement  field. 

We  note  first  that 

7*.  T  =  A(V  •  u)r%  /ir~.  IVu  +  uVi.  (B3) 

The  divergence  of  u  can  be  expanded  as 

V  •  u  =  (rdr  +  r-1  Vj)  *’(T0  +  Vj^)  (B^ 

since 

v  •  (-Tx  VjW)  -  -  T*  (V  X  VjW) 

=  -r*.  (Tx  VdrW) 

=  0 

because  r*  x  V<?rW  must  be  normal  toT.  Equation  B4  can  be  written  as 

V  .  u  =  dtU  +  r-1  V 2 V  +  T«Jr  •  VjK  +  r1  Vj  -Tl/.  (B5) 

The  third  term  vanishes  since  <9r  commutes  with  V*.  and  V^rv  is  normal  to  r.  The  fourth  term  is 
equal  to  (2/r)  U  as  may  be  seen  by  replacing  r_1  Vj  with  V  -  "rc>r,  an  equival&nt  expression.  We 

then  have 
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The  term?*  •  IT*  +  a  Vl  in  eq  B3  is  more  difficult  to  deal  with.  In  terms  of  the  coordinate 
representation  of  *,  (i.e.  uf,  uq,  and  u^), 

T»  IV*  +  *Vl  =  r\ 2dtuf)  +  +  rdf  (^)j  + 

4  'r[nsr<r  V.  *  '*•  [t )J  •  (B7) 

By  inspection  of  eq  5  which  explicitly  gives  the  coordinate  components  of  Vt,  we  see  that  eq  B7 
may  be  rewritten  as 


— *  . 

r 


I  V«  +  ■  Vl  =  7(2dfut)  + 


+ 


+ 


We  note  that  <?r  commutes  with  0  and  <fi.  The  expression  in  square  brackets  represents  the  non- radial 
portion  of  u  and  must  therefore  be  identical  with  Vt  v  -  T  x  Vt  w.  Since  ur  is  identical  with  U, 
we  have 


or 


r  • 


IVu  +  *Vl  =  7te<?rv)  + 


I  Vu  +  u  V  I  =  ~f\2dtu)  + 


Combining  the  above  results,  we  have 


(B8a) 

(B8b) 


(B8c) 
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We  outline,  briefly,  in  this  appendix  the  numerical  techniques  used  in  this  study  to  a)  generate 
the  suite  of  normal  modes  for  a  layered  elastic  sphere  and  b)  evaluate  various  integrals  of  interest 
associated  with  these  modes. 

For  a  given  1  and  w,  the  generation  of  solution  functions  proceeds  exactly  as  outlined  in 
Section  B.  The  matrix  inversion  required  by  eq  45  and  46  was  not  done  explicitly.  We  chose  instead 
to  solve  two  sets  of  simultaneous  linear  equations.  The  Crout  Reduction  (Hildebrand  1956)  was 
found  to  be  particularly  convenient. 

Bessel  functions  were  generated  by  using  Miller’s  well-known  recurrence  algorithm  (Abramowitz 
and  Stegun,  1968).  One  consequence  of  this  technique  is  that  the  accurate  evaluation  of  a  spherical 
Bessel  function  for  many  values  of  its  argument,  as  required,  say,  for  integration,  is  a  time-consuming 
process.  For  such  applications  it  would  perhaps  be  more  efficient  to  numerically  solve  Bessel  s 
equation,  but  computer  memory  limitations  did  not  permit  the  additional  coding  this  required. 

In  practice,  the  program  was  assigned  a  model  and  a  value  of  1  and  proceeded  to  compute  trail 
solutions  for  evenly  spaced  values  of  frequency.  As  the  computation  proceeded,  indicator  variables, 
as  explained  in  Section  B,  were  monitored  for  a  change  of  sign,  which  was  taken  to  indicate  a  zero 
crossing.  When  this  occurred,  an  estimate  was  made  of  the  location  of  the  zero  crossing  and  the 
algorithm  described  below  was  invoked  to  iteratively  improve  the  estir  ue.  In  general,  two  applica¬ 
tions  of  the  following  procedure  sufficed  to  locate  the  eigenfrequency  to  within  one  part  in  10  . 

Gilbert  and  Backus  (1967)  observed  that  Rayleigh’s  principle  could  be  utilized  to  improve 
estimates  of  eigenfrequencies  obtained  by  coarser  methods.  Suppose  that  for  seme  frequency  «, 
near  an  eigenvalue.  «*,  we  have  computed  a  trial  solution  and  find  that  the  stress-free  surface  con¬ 
dition  cannot  be  met.  We  may  apply  Rayleigh's  principle,  or  perturbation  theory,  to  the  solution  we 
have  generated  to  estimate  the  change  in  w  the  elimination  of  surface  stress  would  produce.  The 
first  order  estimate  for  this  change  is  given  by 

_  w.vW  * 1(1  *  (Cn 

5(0  =  (ofTOi)  - — - —  ”  • 

I  ^ 

/  pr2u/f(r)  +  IU  +  l)vf (r)tdr 
o 

We  will  not  derive  eq  Cl  here.  We  then  replace  w  by  w  +  Su  and  repeat  the  process.  We  chose  to 
terminate  the  iteration  when  \8cj/(8o>  +  <u)|  fell  below  10  \ 

The  last  point  we  wish  to  mention  is  the  evaluation  of  integrals  of  the  form 

/  .  J  Z(.)d r  <C2> 

o 
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where  Z  is  some  operator  on  the  solution  functions,  U(r),  etc.  In  general,  we  may  expect  Z  to  vary 
appreciably  over  a  length  L  of  about 


CO 


(C3) 


where  Vp  is  the  local  compressional  velocity.  L  is  simply  the  wavelength  of  a  compression^  wave 
of  angular  frequency  oi.  In  computing  /,  the  program  was  designed  tr  utilize  steps  not  exceeding 
<Lt;  where  t  is  a  small  (~  3  x  10*1)  number  and  Lj  is  the  scale  length  appropriate  to  a  given  shell. 
This  technique  yielded  a  reliably  constant  accuracy  over  many  wide  variations  of  scale  without 
extracting  undue  computing  labor  for  small  values  of  w. 
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APPENDIX  D.  AN  ESTIMATE  OF  TRANSDUCER  AND  SUPPORT 
INFLUENCES  ON  MEASURED  EIGENFREQUENCIES 

Perturbation  due  to  an  anchored  spring 

V/e  suppose  that  one  end  of  a  spring  of  spring  constant  k  is  affixed  to  a  point  P  on  the  surface 
of  the  sphere.  The  spring  is  oriented  radially  and  its  other  end  is  firmly  anchored  at  a  distance 
equal  to  its  rest  length.  When  P  moves  outward  a  distance  u,  it  is  subject  to  a  rorce,  directed 
radially  inward,  of  magnitude  ku.  The  boundary  conditions  for  elastic  motions  of  this  body  car.  be 
expressed  as 

r  •  T  =  -  r  k  (u  •  r  )  8  (r  -  r*p)  on  r  =  rN  .  (°1) 

T  is  the  elastic  stress  tensor,  u  is  displacement,  rp  is  the  location  of  P,  8  is  the  surface  Dirac 
delta  function,  and  rN  is  the  sphere's  outer  surface.  Thus,  for  all  points  on  r  =  rN  other  than  P, 
as  before,  r  •  T  must  vanish. 

We  now  utilize  Rayleigh’s  principle  to  compute,  correct  to  first  order,  the  perturbation  &j“ 
in  an  eigenvalue  u>z  caused  by  the  perturbed  boundary  condition.  The  Fourier  transformed  equation 
of  motion  is 

-  Po>2  u  =  V-  T  (D2> 

(We  have  dropped  sub-  and  superscripts.)  Forming  the  dot  (scalar)  product  of  D2  with  u  gives 

-  pw2 u  •  u  .  (V  •  T)  •  5.  (D3> 

Equation  D3  is  true  everywhere  and  we  may  integrate  both  sides  of  it  over  the  volume  0  <  r  <  rN. 
If  while  doing  so  we  rearrange  the  right-hand  side  somewhat  and  appeal  to  Gauss  theorem,  we  may 
easily  arrive  at 


w2  /  p  u  •  u  dv  /  T  :  V  udv  -  /  r  •  T  ■  udv.  (D4) 

V  v  s 

V  denotes  the  volume  of  integration  and  S  is  its  surface.  The  colon  denotes  a  scalar  tensor  product. 

In  the  absence  of  a  boundary  perturbation  the  surface  integral  vanishes  since  r  ■  T  vanishes  on 
S.  What  remains,  in  this  case,  is  exactly  eq  61  describing  the  equality  of  kinetic  and  potential 
energies  in  free  vibration.  To  determine  the  effect  (to  Hrst  order)  of  the  new  boundary  condition  Dl, 
we  (following  Rayleigh  1945)  replace  «2  by  w2  +  and  evaluate  the  contribution  of  the  surface 
integral.  We  then  subtract  the  unperturbed  energy  terms  (which  are  equal)  leaving 


da>2  f  p  u  ■  udv  k(u  •  r  )2 

v 

where  u  •  r  is  evaluated  at  rp . 


(D5) 
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If  we  adopt  the  normalization  78,  the  volume  integral  becomes  equal  to  unity.  If  rp  is  at  either 
0  or  $  =  n,  and  if  we  replace  /><j2  by  2 wfiw,  we  may  massage  D5  into  the  form 


1 


where  nAj  is  the  appropriate  source  factor,  and  we  consider  only  axisymmetric  (m  =  0)  vibrations. 

A  sphere  agalsst  a  half-space 

In  order  to  estimate  an  appropriate  val'e  for  the  spring  constant  k,  we  consider  the  contact 
problem  of  a  sphere  against  a  half-space.  Let  Eg  and  Vg  be  the  Young's  modulus  and  Poisson's 
ratio  of  the  sphere,  £n  and  Vn  the  corresponding  properties  of  the  half-space,  and  R  the  *adius  of 
the  sphere.  Timoshenko  and  (ioodier  (1951)  present  approximate  results,  useful  when  th>  jontact 
area  is  small,  which  we  quote  here.  The  normal  force  F,  not  stress,  acting  on  the  sphere  causes 
it  to  move  a  distance  d  toward  the  half-space  given  by 


where 


k 


8 


1  -  V2 


(D7) 


and 


k 


n 


'  -  K 


We  suppose  that  the  sphere  is  pressed  against  the  support,  or  transducer,  with  a  force  F  and 
that  motion  results  in  small  variations  aboul  an  equilibrium  position  given  by  D7.  Then  the 
appropriate  value  fur  the  spring  constant  is  given  by 

k  =  WFd)*1  (D8) 


where  the  derivative  is  evaluated  at  F,  the  confining  force.  We  cannot  expect  this  mocsl  to  be  use¬ 
ful  if  any  dimensions  involved  become  significant  compared  to  a  wavelength  of  motion  (about 
2nR/(l  +  Vi).  All  we  really  ask  of  it,  in  any  case,  is  an  order  of  magnitude  for  potential  errors. 
Equations  D7  and  D8  give 


3  16  RF  ft 

2  Ux  ♦  *„)a 


(D9) 


Another  quantity  of  some  interest  is  the  radius  aj  of  the  contact  area,  given  by 
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Table  D1  presents  numerical  results  from  these  expressions  for  a  Lucite  and  a  fused  silica 
sphere.  4  in.  in  diameter,  held  by  its  own  weight  on  a  steel  support  and  a  polyethylene  support. 

The  elastic  constants  were  taken  from  Mason  (1958).  The  lower  portion  of  the  table  gives  the 
effective  spring  constants,  contact  dimension,  and  fractional  perturbations  for  several  modes. 

These  results  suggest  that  soft  materials,  such  as  a  polyethylene,  induce  a  smaller  shift  in 
measured  eigenfrequencies  than  do  hard  ones,  We  may  interpret  this  as  being  due  to  the  low  acoustic 
impedance  of  polyethylene  which  more  nearly  matches  the  free  boundary  condition.  As  impedance 
increases,  the  effective  boundary  condition  in  the  contact  area  approaches  a  rigid  boundary  state, 
and  the  influence  on  eigenfrequency  grows  accordingly.  Equation  D6  requires  Sw/w  to  grow  without 
bound  as  k  goes  to  infinity.  We  must  recall  that  D6  applies  only  to  modes  with  M  =  0.  when  the 
spring  is  affixed  at  one  of  the  poles.  A  spring  affixed  to  a  mode's  node  has  no  disturbing  influence. 

We  also  note  that,  not  surprisingly,  modes  with  small  source  factors  (such  as  jS2)  are  less 
perturbed  than  modes  with  large  source  factors  (such  as  0S2).  Since  the  source  factor  is  directly 
related  to  the  amplitude  of  surface  radial  displacement,  such  a  result  is  to  be  expected.  This  result 
is  not  particularly  helpful  in  choosing  data  since  low-source-factor  modes  are  more  difficult  to 
observe  and  locate  accurately. 
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ratio 

ks  or  kn 

4  In.  sphere 

Lucite 

0.4 

x  10 10 

0.4 

6.68  x  10-" 

8  x  10* 

Fused  silica 

7.3 

x  1010 

0.17 

4.24  x  10'12 

1.38  x  10* 

#347  Stainless 

19.6 

x  10 10 

0.30 

1.48  x  10'12 

Polyethylene 

0.076  x  1010 

0.458 

3.31  x  10-‘° 

Sr 
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pair 

k 

a 

0^2 

0^0 

oS. 

Lucite 

Stainless 

t:  1  x  10* 

8.7  x  10'2 

3.7  x  10'* 

2.5  x  10'5 

2.9  x  10'* 

9.0  x  10'* 

Lucite 

Polyethylene 

2.5  <  10* 

0.16 

1.1  x  10'* 
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Fused  silica 
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0.18 
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1.9  x  10'* 

1.1  x  10'5 

3.9  x  10'7 
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